Monte Carlo Study of Si(lll) homoepitaxy 
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An attempt is made to simulate the homoepitaxial growth of a Si(lll) surface by the kinetic 
Monte Carlo method in which the standard Solid-on-Solid model and the planar model of the 
(7 x 7) surface reconstruction are used in combination. By taking account of surface reconstructions 
as well as atomic deposition and migrations, it is shown that the effect of a coorparative stacking 
transformation is necessary for a layer growth. 
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In contrast to a number of experimental studies on a Si(lll) homoepitaxy [p]-|l6||, there are no theoretical studies 
on it, obviously because of the great complexity of its structure known as the dimer-adatom-stacking-fault (DAS) 
reconstruction jl7]]. Indeed, scanning tunneling microscopy (STM) observations have been the only accessible method 
, to elucidate the growth behavior at the atomic level. 

S In the previous papers, I have proposed the four-state planar model of the Si(lll) surface reconstruction, in which 
the three-point discrete planar rotators (3PDR) were introduced to denote the atoms belonging to the lower half of 
the topmost bilayer (BL). Then, by doing the calculations with the state- flipping dynamics, I have investigated the 
properties of the (7 x 7)-to-"(l x 1)" phase transition fljsfl . Here, the same model is used to examine the growth 
behavior of a Si(lll) homoepitaxy by introducing hopping events as well. This is accomplished by combining the 
3PDR model with the Solid-On-Solid (SOS) model |L9| and let the latter undertake the center force part of the 
interactions. 

In addition to the four states, i.e. unfaulted-stacking (US) states, faulted-stacking (FS) states, and two kinds of 
'-^ \ states which constitute dimer pairs (DA and DB) in Ref. jlq , I am going to introduce a bulk (BK) state, a vacant 
(VC) state, and an intermediate (IM) state which deposited atoms are supposed to take before incorporated into 
bilayer structures. Actually, its importance to dissolve dimer rows and FS halves of the DAS structures were pointed 
out experimentally in relation to the island nucleation and the step flow dynamics [^|8) . Evidently, for these dynamics 
to be realized, the dissociation of FS halves and dimer rows must proceed simultaneously, and thus atoms in the IM 
state must change their states in a coorparative fashion g. To accomplish this, I will choose the atoms belonging to 
the lower halves of the BL's of a Si[lll) surface as the ingredients of the model, and apply the combined model for 
which the algorithm a la Maksym |£0| is employed. Besides, I will assume the FS states to exist only on the topmost 



BL, and prohibit the occurrence of them on underlayers. 

In brief, all terms in the model are described in terms of the 3PDR's states, the basis vectors with i — 0, 1, 2 on 
. the two-dimensional lattice A shown in Figs, [l] (a-b). The first term is, 

oo ; 

H Q =-J ^2 ( 12ce{us,FS} <V(x),<Kx') <W),C ~ <W),us <W'),FS 



x,x'eA 



+ <^(x).US ^(x').BK - <^(x).FS ^(x'),BK j ; (1) 

T3 , 

which generates the US-ordered states, the FS-ordered states, and the bulk structures as shown in Figs. [j] (c-d). In 
this equation, the x and x' summations are taken over all the nearest-neighbor pair sites on A. The symbol y(x) 
denotes the state of a 3PDR at site x on A and ^ denotes the Kronecker's delta. 

The dimer-adatom interaction (DAI) (Refs. pl|-p3[) terms and those connecting them are drawn in Figs. |l| (e-h), 
which are realized respectively by the three-point interactions as 



H\ — — ^l^Z <5 v (x),DA^(x+a;),DB 



xgA i£Z 3 

x [ ^(x-a, + 1 ),US + ^(x-a i+2 ),FS ] i ( 2 ) 

and the terms obtained by interchanging the symbols DA and DB in Eq. (Q). Here, the summation over x is over 
single sites on A, while the equivalence = &j holds for the basis vectors such that i = j (mod 3) in the summation 
over the crystallographic directions (i € Z3). The fact that BK states do not appear in these summations means that 
when an atom sticks onto a US site which is interacting with adjacent dimcrs, the DAI's are diminished or even lost, 
as in reality the dimer contraction M] takes place only when adatoms are associated with dimers. 



1 



The small stacking energy difference between FS and US sites is defined by adding ± J3 to each site, respectively. 

Since I will consider hopping events as well as state-flipping ones, corner holes in the DAS structure must be properly 
dealt with, as opposed to the previous studies in which the center atoms of corner holes remained on a lattice p8| . 
Thus, the interaction term has to be modified to make the central site of a corner hole vacant, which is achieved by 

Hi = ~'^4^ ^ ^(x),DA^(x+a,),DB 

xeAiez 3 ^ 



x 



[^(x-a, t i),VC ^(x-a i+2 ),US + ^(x-a i+ i),FS ^(x-a i+! ),VC 



+ ^(x),DB^(x-ai),DA 

[^(x+a i+ i),VC <^(x+a i+2 ),FS + ^( x+a , +1 ),US ^(x+a, +2 ),VC ] ) • (3) 



X 



The diagrammatic representations of them are shown in Figs. [I] (i~j), in which the cross symbols denote the vacancies. 

Finally, I will introduce an additional term to compensate for unexpected contributions which may arise from the 
SOS model part, because the center force tends to fill in the central site of a corner hole with an unnecessary atom. 
This compensation term is given by, H 5 = 2J 4 J2 xe A Eiez 3 (<W),us + <W),FS ) < W+tu),DA < W-m),DB - 

Thus, the total Hamiltonian is given by the sum of the terms explained above (3PDR part) plus the SOS part 
as H tot = Hgos + kHspdh, where k denotes their relative weight. Correspondingly, the kinetic barrier for the 
growth simulation is defined by — H* t(x) if Wtot(x) < and zero otherwise. Here, 7Y tot (x) is defined from H to t by 
Htot = SxeA^*°*- Then, I will consider as kinetic events the deposition of atoms, a state-flipping event, and a 
hopping of an atom followed by a flipping of its own state. 

By the STM observations, the back-bonding energy of silicon atoms on a Si(lll) surface was estimated to be 
1.020 eV However, since the atoms belonging to the upper half of a BL are also implicitly taken into account 
in the model, it is necessary to include their effect on the back-bonding energy of the SOS part J su b as well to give 
it a slightly larger value. On the other hand, Ichimiya and co-workers estimated the activation energies of the island 
decay and hole-filling phenomena to be around 1.5 eV and 1.3 eV, respectively [[ll 12 . As opposed to the former 



case, these values should be larger than J sv b because other interaction effects are included in their measurements. 
Therefore, I will take J su b = 1.20 eV, which is near the average value of the reported ones. 

At high temperatures, surface structure formation will be dominated by the SOS part, and thus the coupling 
constant for the lateral part Ji a t will give the principal contribution to the surface melting temperature ss 1400K. 
Therefore, it is natural to set Ji a t = 0.12 eV. 

As for the 3PDR part, the ratios Ji/Jo (i = 1, . . . , 4) of the coefficients of the interaction terms Hi appropriate in 
producing the (7 x 7) DAS structure were already known |T^] . To be precise, J4 has not yet been estimated, because the 
Hi term in Eq. (||) takes the different form than in Ref. JTqj. However, this difference arises merely due to the absence 
of hopping processes in the former work, and the same value for J 4 should be applicable also to the present case. The 
ratios of the parameters thus used in the calculations are given by, J1/J0 = 4.175, J2/J0 = 2.450, J3/J0 = 0.100, 
and J4/J0 = 0.500. 

Now, the only remaining unfixed parameter is k. In the previous study Jig ], it was also shown that the transition 
temperature T t of the (7 x 7) reconstruction was given approximately by Jo, so that kJq » fcsTj, where ks is the 
Boltzmann's constant. Then, since T t has been reported to be about HOOK p4[ ), I will set kJq — ksT t = 9.5 x 10~ 2 eV. 

It has been shown experimentally that the complexity and the stability of the (7 x 7) structure obstruct the growth. 
Moreover, when atoms are deposited, the kinetic process has been shown to proceed by first forming the clusters of 
IM states and later by their transformation into US states To accomplish this, I will introduce the interaction 

energies of this state, i.e. I will denote the lateral interaction between a pair of IM states by E* at , and its back-bonding 
energy on a US site by E^g and that on a FS site by Ep S , respectively. 

Suppose that an atom is going to stick onto a site x with its three substrate atoms in a plaquette at x, x + a.2 , and 
x — a.±. Then, the new state of an atom after a deposition or a hopping is determined by the following rules: (1) if all 
three substrate atoms are in the US states, then set it to be in the US state if any one of the adjacent atoms in the 
same plane is in the US or in the bulk state, and if not, set it to be in the IM state if any one of the adjacent atoms 
in the same plane is in the IM state or it is isolated, (2) if some of the substrate atoms are either in the US states or 
in the FS states but both are not simultaneously present, then set it to be in the IM state if any one of the adjacent 
atoms in the same plane is in the IM state, (3) if all three substrate atoms are in the FS states, set it to be in the IM 
state, and (4)If none of them is met, an atom is prohibited to stick there. Here, due to the bilayer character of the 
3PDR's, I will assume no bonding is available for an IM-state atom to prevent further sticking of an atom on it. If 
the size of an IM-state cluster exceeds a critical size S cr , a stacking transformation takes place and all atoms in the 
cluster is set to be in the US states. 
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The length of the random hopping process is set to be 50, so that the effective radius of the searching area 
approximately becomes equal to the linear dimension of the (7 x 7) unit cell. Then, by carrying out the calculations, I 
found it appropriate to choose E* at = 1.35 eV, E^ s = 1.20 eV and E FS = 1.30 eV to produce an oscillatory behavior 
for the step density (SD) when the atomic flux and the substrate temperature are chosen to be 0.1 monolayer coverage 
per second and 650 K, respectively. In Fig. |^, (1 - SD) is plotted as a function of a coverage for which the hopping 
length is 50, and S cr = 17 is used. The lattice sizes are 35 x 35 and 49 x 49, where the snapshots of the latter are 
shown in Fig. [|. 

Interestingly, the oscillatory behavior of the SD was seen only when kJq is very close to 9.5 x 10 -2 eV. Unfortunately, 
however, this does not look like the ordinary BL-by-BL growth behavior observed by reflection high-energy electron 
diffraction (RHEED) experiments and by STM observations J|,|l5],[l6| . Instead, it looks similar to the RHEED intensity 
oscillation obtained for a Ge/Si(100) heteroepitaxy p5| . This means that atop sites of a growing surface is favored 
for further growth. At the same time, this island-growth behavior appeared due partly to the limitation of the 3PDR 
model, with which only the bilayer character of the growth can be realized, and hence this discrepancy implies the 
true IM state does not possess this character. In other words, with the present model, there is no guarantee that the 
stacking transformation into the US ones is irreversible, and this property seems to play a very crucial role in the 
surface growth. Actually, in the present model, it is not possible to preclude an event in which an atom in a US state 
to hop around and stick to an IM island and change its state into an IM one. 

In summary, by the combined use of the 3PDR model |l8| and the SOS model p9| , the kinetic growth simulation of 
Si(lll) homoepitaxy is performed. It is found that atop sites of a growing surface is favored for further growth. Also, 
it is found that for the bilayer step density to show the oscillatory behavior, it is necessary to choose the back-bonding 
energies so that the inequality J su i, = E^ s < Ep S holds. However, this result is evidently related to the bilayer 
character of the 3PDR model. More importantly, the irreversibility of the transformation from an IM state to a US 
one is crucial in obtaining the BL-by-BL growth mode. 
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FIG. 1. Interaction terms of the Hamiltonian are displayed diagrammatically. (a) Allowed four states of the 3PDR's, 
US,FS,DA, and DB. (b) Three basis vectors used in the equations, (c-d) Locally US and FS ordered states, (e-f) DAI terms, 
(g-h) DAI-connecting terms, (i-j) Corner-hole stabilizing terms. 

FIG. 2. (1 - SD) plotted against the coverage. The lattice sizes are 35 x 35 and 49 x 49, the hopping length is 50, and the 
critical cluster size is 17. 

FIG. 3. Snapshots of simulation on Si(lll) homoepitaxy. The lattice size is 49 x 49. Atoms at higher layers are depicted by 
darker and larger symbols, where 3PDR's denote US,FS,DA, and DB states, and shaded discs denote IM states, (a) 0.40 BL. 
(b) 0.70 BL. (c) 1.00 BL. (d) 1.50 BL. 
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M.Itoh Figure 3(c) 




M.Itoh Figure 3(d) 



